function [nbrmat] = findnbrs(imat,nbrs)

[m,n] = size(imat);
nbrmat = zeros(m*n, nbrs);

for j=1:n
    for i=1:m
        % above
        if iflegal(i-1,j,m,n)
            nbrmat(sub2ind([m,n],i,j),1) = sub2ind([m,n],i-1,j);
        else
            nbrmat(sub2ind([m,n],i,j),1) = 0;
        end
        %left
        if iflegal(i,j-1,m,n)
            nbrmat(sub2ind([m,n],i,j),2) = sub2ind([m,n],i,j-1);
        else
            nbrmat(sub2ind([m,n],i,j),2) = 0;
        end
        % below
        if iflegal(i+1,j,m,n)
            nbrmat(sub2ind([m,n],i,j),3) = sub2ind([m,n],i+1,j);
        else
            nbrmat(sub2ind([m,n],i,j),3) = 0;
        end
        % right
        if iflegal(i,j+1,m,n)
            nbrmat(sub2ind([m,n],i,j),4) = sub2ind([m,n],i,j+1);
        else
            nbrmat(sub2ind([m,n],i,j),4) = 0;
        end
    end
end

nbrmat = nbrmat';

function ifok = iflegal(i, j, nr, nc)
if i>=1 && j>=1 && i<=nr && j<=nc
    ifok = 1;
else
    ifok = 0;
end